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Introduction 

In  this  report  we  consider  the  general  problem  of  contact  and  impact 
between  two  bodies.  The  report  is  divided  into  three  basic  parts.  These 
parts  describe:  (I)  The  general  theory  of  contact-impact  problems, 

(II)  A numerical  scheme  for  the  analysis  of  contact- impact  problems,  and 

(III)  The  description  of  computer  program  FEAP  74  for  the  solution  of 
contact- impact  problems.  In  an  appendix  we  include  the  program  subroutines 
and  general  input  description  for  FEAP  74. 

In  Sections  1 to  6,  Part  I,  we  deal  with  spatial  aspects  of  the  theory 
and  in  Section  7,  Part  I,  we  deal  with  temporal  aspects.  This  splitting 
of  the  theory  is  motivated  by  the  way  we  intend  to  numerically  solve  the 
equations,  i.e. , the  finite  element  method  spatially  and  a finite  difference 
method  temporally. 

Part  II  considers  a numerical  implementation  of  the  theory  given  in 
Part  I.  Section  9 deals  with  spatial  notions  of  the  numerical  problem 
and  Section  10  the  temporal.  The  solution  scheme  for  the  resulting 
algebraic  problem  is  discussed  in  Sections  11  and  12. 

The  computer  program  FEAP  74  was  modified  to  incorporate  the  numerical 
contact- impact  model.  The  program  modifications  and  capabilities  together 
with  two  numerical  examples  are  contained  in  Part  III. 

Finally,  in  the  appendix  we  give  listings  for  the  contact  subroutines 
together  with  the  data  input  instructions. 
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PART  I 

VARIATIONAL  FORMULATION  OF  CONTACT- IMPACT 
PROBLEMS  IN  CONTINUUM  MECHANICS 


1 . Prel iminaries 

Our  conventions  on  indices  are  as  follows: 

Superscripts  indicate  to  which  body  an  entity  pertains.  Sumnation  is 

to  take  place  only  when  explicitly  indicated. 

Latin  subscripts  range  over  1,2,3,  while  Greek  subscripts  range  over 

1,2.  The  summation  convention  is  assumed  to  hold  for  both. 

A body  i s a nice  connected  region  of  [Rswith  a piecewise  smooth 

★ 

boundary  d<©  . A contact  problem  is  a boundary  value  problem,  or  an 
initial -boundary  value  problem,  in  which  two  bodies,  ® and  <8 , interact 
according  to  the  principles  of  mechanics.  Thus  the  primary  kinematic 

Lilt 

axiom  of  a contact  problem  is  that  configurations  .6'  and  (r , of  ft  and  ft  , 
respectively,  do  not  penetrate  each  other,  i.e., 

CJr‘)'n  V - 

_jr*-  r\  c -V)"  « 

where  denotes  the  interior  of  -tT  , «*  = 1,2. 

On  the  other  hand  the  unique  conditon  which  characterizes  contact 
problems  is  that  material  points  on  the  boundaries  of  •'->  and  ~ may 

1 t 

coalesce  during  the  motion  of  the  bodies.  Thus  we  say  j and  'ft  are  in 
contact  If  dJr1  & , and  we  define  the  contact  surface 

e by 

_ 

It  is  usual  for  the  term  contact  to  have  a static  connotation  while  the 
term  impact  has  a dynamic  connotation.  We  shall  use  contact  in  the  general 
sense  to  include  static  as  well  as  dynamic  phenomena. 
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te.  - n dJj1  . (2) 

If  <£>  and  ® are  never  in  contact  then  t = <&  for  all  configura- 
tions Jr’andJ^,  and  in  this  case  an  initial -boundary  value  problem  for 
t^and  ^ reduces  to  one  in  w.i:>  h & and  ^ may  be  treated  separately. 
Thus  a non-trivial  contact  problem  is  one  in  which  for  at 

least  one  instant  during  the  motion  of  S^and  The  picture  (Fig.  1) 
illustrates  these  notions. 

Equation  (1)  implies  that  e is  a material  surface  with  respect  to 
both  bodies,  i.e.,  one  which  is  not  crossed  by  material  particles.  From 
this  we  may  deduce  the  interface  conditions  on  te  . 

Let  x be  a persistent  point  of  fe.  (one  at  which  joining  or  re- 
leasing of  the  bodies  is  not  instantaneously  occurring)  and  v be  the 
velocity  of  x ( y = x ).  Note  that  only  the  normal  part  of  y is 
independent  of  the  parametrization  of  fc  . Let  y1  and  y*  be  the 
velocities  of  the  material  particles  located  at  the  points  x1  and  **, 
contained  in  db1  and  dir1  , respectively,  such  that  5 =;‘  = x* 
at  the  present  instant.  Then  since  fe  is  material  and  x is  persistant 

V-  n « vl-  n « v*  n (3) 

where  2 is  a unit  normal  vector  to  fc.  at  x . From  this  it  follows  that 
a necessary  condition  for  momentum  to  be  balanced  at  x is  that 

( tl  - t‘)  n - o , (4) 

where  t*  is  the  Cauchy  traction  vector  with  respect  to  tJbr* 

In  addition  we  assume  that  no  tensile  tractions  can  occur  on  e.  , 
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i*  *2  *=  ° , (5) 

where  n“  is  the  outward  unit  normal  vector  to  dir*  This  condition 
excludes  the  possibility  of  the  two  bodies  being  glued  together.  Condi- 
tions (1-5)  characterize  our  notion  of  a contact  problem. 

Note  that  thus  far  we  have  said  nothing  about  the  tangential  parts 
of  v ' and  t**  . These  remaining  conditions  are  determined  by  the 
frictional  nature  of  the  contact.  We  shall  study  two  simple  cases. 

Cas(M:  If  we  assume  that  points,  once  in  contact,  move  with  fc 

until  released,  we  have  that 

, (6) 

and  therefore 

t1  ♦ t‘  = o . (7) 

For  this  model  we  say  that  a no-slip,  or  perfect  friction,  condition  is 
achieved  on  fc  . Thus  condition  (5)  and  equations  (6)  and  (7)  are  the 
interface  conditions  for  this  case. 

Case  II:  We  may  create  the  interface  conditions  for  a frictionless, 

sliding  contact  by  asserting  that  the  tangential  part  of  each  t"  is 
identically  zero, 

f - (t  n-)n-  . c (d) 


Eq.  (8),  along  with  (3-5),  are  the  interface  conditions  for  this  case. 
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2.  Variational  Theorems 

We  will  formulate  a variational  theorem  for  the  contact  problem  of 
finite  elastodynamics . We  point  out,  however,  that  our  treatment  is 
entirely  general  and  could  be  used  in  conjunction  with  any  field  theory, 
as  the  only  unique  feature  of  the  formulation  involves  the  handling  of 
interface  conditions.  At  the  same  time  finite  elastodynamics,  though 
lending  itself  to  a clean  and  simple  variational  statement,  is  a case  of 
wide  practical  interest. 

We  shall  first  obtain  a variational  theorem  for  the  usual  initial- 
boundary value  problem  of  finite  elastodynamics  by  a trivial  generaliza- 
tion of  some  work  done  by  S.  Nemat-Nasser  [1]. 

For  notational  simplicity  let  ci.  denote  dth,  and  let  cut  and  dfi  denote 
area  and  volume  forms  for.^andCL,  respectively.  Let  CJ^c  GL  be  that 
part  of  Ct  where  surface  tractions  are  prescribed,  and  denote  by  T the 
Piola  - Kirchhoff  traction  vector  representing  these  prescribed  tractions. 
Call  p.  the  density  of '13  in  the  initial  configuration,  F the  extrinsic 
body  force  vector  and  let  * - *,(*)  represent  the  position  at 
time  t of  the  material  particle  located  at  X in  the  initial  configuration. 
For  convenience  we  take  to  be  the  initial  configuration.  We  denote  by 

(* X the  deformation  gradients  and  by  <£(<?¥ /<3X  ) the  strain 
energy  density.  Then  if  x satisfies  the  kinematic  boundary  conditions 

* = (9) 


on  CL  , where 


a.  U aT  = a , 

naT  = 0 , 


the  functional  IE  defined  by 
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HCs)  = 5>  f ~ 

° <Q 

- p.  F x ) <±&  - ^ ■*  • T d.CL  | d.t  , 


(10) 


is  stationary,  i.e.,  its  first  variation  vanishes 

o = ) * Sf  H p.  (s  - f ) - D,v  ? >**«*a 

• e 

■«-  ^(T-T)ix  cLGt  | dLt  , 
dr 


(ID 


subject  to  the  constraint  on  variations  Sx  - W,  = o > (12) 
if  and  only  if  the  equations  of  motion  and  traction  boundary  conditions 
are  satisfied 


f-  ( 2 “ f ) = 

DIV  P , 

in  6 j 

(13) 

T *■ 

T , 

on  dT  f 

(14) 

where  P~  d<p/d(d*/dX  ) is  the  first  Piola  - Kirchhoff  stress  tensor, 
T = N P is  the  Piola  - Kirchhoff  traction  vector,  and  M is  the  out- 
ward unit  normal  vector  tou‘.  The  solution  to  the  initial-boundary  value 
problem  must  also  satisfy  the  given  initial  conditions 


x » x, 

Y = Y- 


in  43 


(15) 


To  interpret  this  variational  theorem  for  two  (non-interacting) 


e - 

a - 

a1  ^ 

a*  , 

etc 

bodies  set 


8 


and  write 

Xfe)  = +■  ffV**)  . 


The  next  step  is  to  add  to  I terms  manifesting  the  interface 
conditions  on  fe.  and  to  stipulate  the  constraints  under  which  the  vanishing 
of  the  first  variation  of  the  appended  functional  corresponds  to  a 
solution  of  the  contact  problem.  To  do  this  we  must  consider  further  the 
kinematics  and  geometry  of  "e.  . 

Define  two  piecewise  smooth,  invertible  maps  by  the 

condition 


£ ^ a" 


(16) 


where  each  a*  identifies  points  on  the  boundary  of  the  initial  con- 
figuration £0“  which  map  into  the  contact  surface  'r.  at  each  instant  of 
time.  If  , then  X ^(^  1 ) Y*  ) and  X * = (£  \)  Y*  ) are  the 

positions  of  particles  in  GL  a ndt<.\  respectively,  which  have  coalesced 
at  *€  fc  . It  is  clear  what  the  ^"*'s  really  are,  viz.,  if  « 

for  all  x%  O*4  , represents  the  motion  of  body  £6“  from  the  original 
configuration  <fT  to  the  present  one  A' , then  is  the  restriction  of 
> to  £ , 

*'«■)  - *"(*')  , (,7) 

for  each  x"e  , *=1,2.  . For  the  time  being  we  consider  the  -jC'  's 

as  maps  defined  independently  of  the  v*4  's  and  consider  (17)  a constraint 
on  possible  motions. 

We  are  interested  in  to  what  extent  the  relation 
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(18) 


is  smooth  in  time  and  analogously  under  what  circumstances  the  variations 
of  the  •/.*  's  are  equal.  In  general  the  -/■*'  's  will  not  even  be  contin- 
uous in  time  since  contact  surfaces  can  be  instantaneously  created  «r 
destroyed.  If  we  eliminate  such  exceptional  instants  and  consider  only 
persistent  points,  the  bodies  still  may  slide  with  respect  to  each  other, 
as  depicted  in  Fig.  2.  Thus  tangential  velocities  are  seen  to  be  unequal 
in  general.  However,  when  * is  persistent,  the  impenetrability  condition 
(1)  forces  the  normal  velocity  components  to  be  equal,  and  concomitantly 
the  normal  components  of  variations  of  the  ^ 's  are  also  equal 


(19) 


For  sliding  contact  (Case  II),  Eq.  (19)  characterizes  the  constraint  on 
variations  of  the  's  eqivalent  to  the  velocity  constraint  (3). 

For  no-ilip  contact  (Case  I), 


s^1  - , (20) 

is  easily  seen  to  be  the  condition  on  variations  equivalent  to  Eq.  (6). 
We  shall  see  that  Eqs.  (19)  and  (20)  lead  to  the  proper  interface  condi- 
tions in  the  variational  theorems. 

Introduce  vector  valued  Lagrange  multipliers  , and  add 
2 * 

\ \ T-.  d.£-d.t  , (21) 

o £«  - 

to  the  functional  I (Eq.  10).  Note  that  when  0 , 


« OL~  u oT_  u 

■x  T > 


a 
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Figure  2 
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and  assume  for  consistency's  sake  that 

T a O 


(22) 


This  condition  will  preclude  the  ambiguous  circumstance  of  non-zero 
tractions  being  specified  on  the  contact  area.  Upon  taking  variations  of 
$-  = ![-*- X we  get  Eqs . (13),  (14)  and. 


o = 


ai  n } C x“-  -JT  ) 


(23) 


(T~-  T~)  - T1"*  jcLt'dit 


* transversal  i ty  condition. 

The  transversal ity  condition  is  the  classical  terminology  for  variations 
associated  with  the  domain  (Z 

The  first  summand  of  (23)  gives  us  (17)  which  insures  that  the  's 
map  into  fc.  properly.  The  second  suirmand  identifies  Tr  as  the  Piola  - 
Kirchhoff  traction  vector  T "‘on  C" . Let  us  investigate  the  third 
summand. 

Consider  first  Case  I and  define 


= S 2^  , ocu.i  j (24) 

which  makes  sense  because  of  Eq.  (20).  This  condition  is  equivalent  to 
insisting 

> ■*-  > 


thus  the  first  sunmand  of  (23)  also  implies  (6)  holds  whenever  we  have  a 
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persistent  point.  Let  denote  the  Jacobian  determinant  associated 
with  ^ , 

de  ~ (25) 

Notice  then  that  since  is  the  Piola  - Kirchhoff  traction  vector, 

c*  is  the  corresponding  Cauchy  traction  vector.  With  these  we 
have  for  the  third  summand, 

= ==<<,  iaT s'lC- - \ s^-Ccv^r'.a//)^  )ctt  , (26) 

-_i  <r  * 

which  in  words  means  the  Cauchy  traction  vectors  are  in  equilibrium.  Thus 
the  momentum  balance,  Eq.  (7),  is  satisfied  on  >e  . 

In  Case  II  we  only  have  that  (19)  holds,  so  define 


S^C")  — n , i,2  . (27) 

This  requirement  also  insures  that, 

v* 

* 0 - * ‘ C , 

thus  the  first  summand  of  (23)  implies  (3).  For  this  case  the  third 
suitmand  takes  the  form, 

C • (28) 

W(">a) 

•**!  £•* 

The  integral  over  fe  gives  us  Eq.  (4).  The  significance  of  the  second 
integral  hinges  on  the  observation  that  S^cnjn  ) is  a tangent 

vector  to  for  each  . Thus  the  tangential  part  of  each  T*  is 
identically  zero,  which  is  equivalent  to  the  shear  free  condition,  Eq.  (8), 
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which  we  require  for  Case  II. 

A standard  calculation  enables  us  to  write  the  transversal ity 
condition  as, 

O -ks  CSX"  'T'xf  ■(«--  ■4-'))  teat-)  , (29) 

•**  1 w ^ 

at 

where  the  transversal  T*  is  a unit  vector  field  tangent  to  , and 
perpendicular  ard  pointing  outward  with  respect  to  dt"*  , Fig.  3.  Thus 
(29)  implies  that 

C ) = o on  , -K-i.i  ^ 30) 

Assuming  continuity  of  the  integrands  of  (21)  on  the  closure  of  C , 
condition  (30)  is  already  implied  by  the  first  summand  of  (23).  This 
assumption  precludes  taking  the  form  of  a & -distribution  on  . 

Although  this  assumption  is  warranted  here  it  may  not  be  true  when 
one  employs  certain  approximate  theories  in  mechanics.  For  instance 
consider  the  case  where  a Bernoulli -Euler  beam  is  uniformly  loaded  and 
sits  on  u rigid  parabolic  surface  (Fig.  4).  At  the  contact  points  a,  a', 
concentrated  reactions  must  exist  to  balance  shear  forces.  This  example 
is  actually  from  a completely  different  class  of  contact  problems  in  that 
contact  is  made  along  a part  of  the  interior  rather  than  the  boundary. 

Such  problems  as  the  contact  of  plates  and  shells  also  fall  into  this 
class.  We  could  summarize  such  situations  by  the  description  — 
m-dimensional  contact  of  m-dimensional  bodies,  e.g.,  for  the  beam  m=l, 
and  for  plates  and  shells  m«2.  The  case  under  investigation  in  this  paper 
(m«3)  is  an  example  of  the  (m-1 )-dimensional  contact  of  m-dimensional 


bodies . 
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Figure  4 
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It  is  good  to  keep  in  mind  cases  such  as  that  illustrated  in  Fig.  4 
when  considering  specific  boundary  value  problems. 

A further  point  worth  mentioning  here  is  that  the  transversal ity 
condi t an  will  in  general  be  an  Independent  one  in  a numerical  algorithm. 
For  example,  if  the  fields  in  the  integrand  of  (21)  are  approximated  by 
a family  of  trial  functions,  Eq.  (23)  only  implies  that  some  weighted 
integrals  over  the  (2. " 's  vanish.  The  condition  (29)  requires  that 
weiqhted  integrals  over  the  <bC~  ' s also  vanish. 

We  now  summarize  our  results  in  the  following  theorems: 

Theorem  I:  Let  (1),  (2),  (5),  (9),  (12),  (15),  and  (20)  hold.  Then 

x is  a solution  to  the  no-slip  contact  problem  (Case  I),  that  is, 

(6),  (7),  (13),  (14)  and  (17)  also  hold  if,  and  only  If,  S^ofor 
arbitrary  variations  of  x"\  and  *«i,z  . 

Theorem  II:  Let  (1),  (2),  (5),  (9),  (12),  (15),  and  (19)  hold.  Then 

x is  a solution  to  the  sliding  contact  problem  (Case  II),  that  is,  (3), 
(4),  (8),  (13),  (14)  and  (17)  also  hold  if,and  only  if,  = o for 
arbitrary  variations  of  v-,  “ and  'C~  . 
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3.  Consideration  of  Theorems  I and  II  as  Computational  Tools 

Theorems  I and  II  may  be  employed  to  generate  numerical  algorithms 
for  the  solution  of  contact  problems.  The  basic  idea  is  to  represent  -x- , 
a-- and  <T  as  the  product  of  known  functions  on  IR3  with  unknown 
parameters  depending  on  time.  Then  Theorems  I and  II  provide  us  with  a 
method  for  generating  an  approximate  system  of  equations  (e.g.,  by  the 
classical  Ritz-Galerkin  technique)  in  terms  of  these  unknown  parameters, 
which  then  can  be  solved  incrementally  and/or  iteratively,  subject  to  the 
side  conditions  of  tie  theorems.  The  constraints  (1)  and  (5)  will  both 
take  the  form  of  inequalities  in  actual  computations,  this  the  ideas  of 
optimization  theory  will  probably  be  useful  in  the  actual  construction  of 
a numerical  algorithm. 

The  finite  element  method  is  a powerful  technique  for  obtaining  a 
system  of  approximate  equations,  and  it  is  of  interest  to  find  out  how 
amenable  are  Theorems  I and  II  to  a finite  element  formulation. 
Unfortunately  the  term  3^  would  result  in  a terrible  mess  if  the  integrand 
was  represented  by  typical  finite  element  functions.  This  is  because  the 

to 

boundaries  of  the  C 's  are  unknown  and  thus  a parametric  integration 
would  bury  the  defining  parameters  of  the  C ' s in  the  arguments  of 
Heaviside  functions  representing  the  supports  of  the  elements.  Note  that 
a classical  Ritz-Galerkin  approximation  would  not  be  subject  to  this  pit- 
fall,  since  the  associated  trial  functions  could  be  chosen  to  be  real 
analytic  and  thus  easily  integrated  parametrically  to  a relatively  simple 
form.  However,  such  a formulation  is  restricted  to  a geometrically 
simpler  class  of  problems.  Thus  it  is  desirable  to  seek  a generalization 
that  will  lend  itself  cleanly  to  a finite  element  formulation. 
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4.  Variational  Theorems  Without  Transversal ity  Conditions 
Let  C be  a fixed  part  of  such  that 

e-  => 


(31) 


and 


T = o 


(32) 


Define  a scalar  valued  function  -t?"  on  such  that 


V(x-)  = ° if  x“t 


Let  e j>  r.  , and  define  the  maps  ^ by  the  condition 


L2Sr  ) 


(33) 


where,  as  before,  £ represents  / on  C ; but  on  C ~~  C we  place 

, •»  *■< 

no  physical  interpretation  on  . Thus  on  C we  will  always  have  that, 

c 2^")  - s , (34) 

since  x**-  ^ on  C and  >2*°  on  the  relative  complement 

Introduce  vector  valued  Lagrange  multipliers  and  let^t=I*  ni 

where 


TTt  ^ i \ \ <r~  .jlC' 

o ^ 


i t 


(35) 


We  require  that  the  variations  of  ■)/.  * satisfy  the  same  conditions  as 
before,  but  now  for  all  (Z*  : 


Case 


f . 


1 


Case  II:  ^/(n)  n 


on  C. 


(36) 
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where  n is  a unit  normal  vector  to  e . Computinq  the  first  variation 
of<£  we  have  the  usual  conditions  emanating  from  31  and 


= - =£  S S { i0-**  CVc  ^m)) 

o jg  — L 


(37) 


Sx“  CrzV  - T - ) 


-*i/-(Vr-)  )*£■ 


cLt 


The  first  summand  gives  us  (34)  and  we  define 

t.'“-  ( x-t£“  = «‘(x“)-;A**)  } . (38) 

The  third  sumnand  defines  “*?**  or  “ as  the  Piola  - Kirchhoff  traction 
vector.  Note  that  this  insures  that  on  C^~  C*  since  >Z^  ° there. 

The  fourth  sumnand  gives  us  the  appropriate  Cauchy  traction  condition 
across  e for  each  case  of  (36).  The  second  summand  is  identically 
satisfied  on  C since  . On  t — C-  it  tells  us  that 

is  orthogonal  to  , but  this  is  of  no  physical  interest 

Thus  we  can  state  the  following  theorems: 

Theorem  1 1 : Let  (1),  (2),  (5),  (9),  (12),  (15),  (31),  (32)  and  (36)] 

hold.  Then  > is  a solution  to  the  no-slip  contact  problem  (Case  I),  that 

is,  (6),  (7),  (13),  (14)  and  (17)  also  hold  where  t is  defined  by  (38), 

if  for  arbitrary  variations  cf  and.  cr  “ , - - l , z 

Theorem  IT:  Let  (1),  (2),  (5),  (9),  (12),  (15),  (31),  (32)  and 

( 36 ) ^ hold.  Then  ^ is  a solution  to  the  sliding  contact  problem  (Case 

II),  that  is,  (3),  (4),  (8),  (13),  (14)  and  (17)  hold  where  C"*  is 

defined  by  (38),  if  ° for  arbitrary  variations  of  ^"*,“7 

<r  ~ * i . z 
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Thp  important  feature  of  these  theorems  is  that  the  regions  ^ 
are  fixed.  Thus  transver sal ity  conditions  are  absent,  and  the  theorems 
may  be  applied  to  finite  element  formulations.  In  fact  one  would  naturally 
take  C to  be  a union  of  elements  in  , large  enough  to  contain 

C.  throughout  the  motion. 

Thus  far  our  considerations  have  been  quite  general  and,  in  fact, 
i, tore  qeneral  than  would  be  required  for  the  solution  of  particular  classes 
of  contact  problems.  In  the  next  section  we  illustrate  the  many 
simplifications  which  can  be  made  in  the  application  of  the  preceeding 
theorems  to  a class  of  problems  of  wide  practical  interest. 
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5.  Hertzian  Contact  Problems 

We  wish  to  characterize  contact  problems  in  which  the  contact  surface 
is  approximately  planar  and  the  bodies  have  undergone  small  deformations 
in  the  neighborhood  of  the  contact  surface. 

Assume  the  following: 

(1)  2^frt,e,  58  Ss  on  fe  , where  the  indicate  components 

. . 3 3 

with  respect  to  the  standard  basis  [ e.  j for  |R  , 

(see  fig.  5). 

(2)  1 . -*-1,2.  , thus  X on  . 

Assumptions  (1)  and  (2)  together  imply  that, 

«*  _ * 

' G ~ X C T 3 , 


and  that, 

(t“,  C.°  ) ~ £**  Ct‘n)n  = I--  (I  n)n  « C.  \ > Tt  , O ) . 

(3)  Material  points  which  eventually  contact  have,  to  the  first 
order,  the  same  initial  coordinates  and  hx  . Explicitly 
we  manifest  this  idea  by  requiring  that  the  's  satisfy 

& ^ X *(«*.«,))  . (39) 

This  i„  depicted  in  Fig.  6.  Since  X*  are  given  functions  which 
define  the  surfaces  t , it  follows  from  (39)  that, 

= i/:1  . 

We  term  problems  for  which  these  assumptions  hold  Hertzian,  since 


these  assumptions  are  implicit  in  Hertz'  classical  theory  [2]  (see 


Figure  5 
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Goldsmith  [3]  for  an  excellent  exposition  of  this  work  and  also  many 
applications  of  hertz'  theory  to  impact  problems).  It  should  be  pointed 
out  that  the  formulation  we  are  about  to  give  is  still  considerably  more 
general  than  thnse  to  which  Hertz*  theory  applies. 

We  now  sh  v how  these  assumptions  allow  us  to  make  simplifications 
in  the  precedi n„  theorems. 

Theorems  I and  II: 

Due  to  assumption  (3)  the  term  can  be  replaced  by  an  integral  over 

a region  in  the  a^z,.  -plane.  This  region,  say  c , is  the  projection 
of  £ onto  the  z^.z^  -plane,  and  due  to  assumption  (2)  it  coincides,  to 
the  first  order,  with  the  projections  of  the  's.  Thus  3<lcan  be 

written 

(40) 

e,  c 

Since,  for  Case  I,  we  know  that  the  momentum  balance  on'e.  requires  that 


we  may  make  use  of  this  relation  immediately.  Thus  define 


and  substitute  into  (40).  Employing  (39),  the  integrand  simplifies  to 

T'  (**-  (41 ) 


The  analog  of  (23)  becomes 
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© = S S { sr  c*2-  x1)  - 

o C.  I 

«-  - ^5*  ) J dcdt 

(42) 

-+  transversality  condition  . 


Thus  the  same  conclusions  of  Theorem  I can  be  drawn.  However,  fr«m  a 
numerical  standpoint  things  are  considerably  different.  First  of  all, 
since  the  jt*  's  are  absent  in  this  formulation,  we  do  not  get  a uniquely 
defined  b ; x1  and  x*  will  not  in  general  be  the  same  pointwise.  If 
the  graph  of  b is  important  it  could  be  constructed  by  averaging  * 1 
and  x *■  , which,  if  the  solution  is  any  good,  should  be  reasonably  close 
pointwise.  On  the  other  hand,  the  ;/"*  's  being  absent  engenders  a con- 
siderable saving  in  the  number  of  equations  to  be  solved  and  in  their 
complexity. 

The  analogous  case  for  Theorem  II  is  constructed  simply  by  setting 
Then  the  integrand  of  yC  becomes 


T (x 


1 

5 


and  (42)  reduces  to 

° - S S [ stc*,’  -*.1)  ♦ 

o c 

- ix,1  T.1  4 ix.1  Tj1  j dcdt 


4 transversal) ty  condition 


(43) 


(44) 


Hence  the  conclusions  of  Theorem  II  hold. 
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Thus  in  the  case  of  Hertzian  contact  we  can  add  the  simplifications 
manifested  in  (41)  and  (43)  to  the  conditions  of  Theorems  I and  II, 
respectively,  and  still  gamer  the  same  conclusions. 

Theorems  1 1 and  II 1 : 

For  these  cases  *771  can  be  written  as  an  integral  over  c , the 
projection  of  te  : 

^ S S die  dLt  . 

Due  to  the  present  geometric  situation,  ic  is  appropriate  to  take 


and  thus  define 


- ~>Z  > -*-1,2. 


Analagous  to  the  considerations  for  Theorems  I and  II,  the  momentum 
balance  across  € motivates  the  simplification 


*Mi  i.  I 

<t  — <t  * _ <r 


With  these  and  (39),  the  integrand  of  771,  can  be  written 

£ • C x - i ) ■ 

* 

A further  simplification  can  be  made  by  setting 


This  eliminates  one  unknown  function  and,  as  we  shall  see,  has  the  effect 

of  satisfying  (5)  naturally.  Thus  the  integrand  of  "777  becomes 
- ’ — 1 

This  is  a standard  ploy  of  optimization  theory,  see  p.  82,  [4]. 


27 


OU  ~ **  ) - <1^2  )iC 


- V, 


(45) 


and  the  analog  of  (23)  is 

° = S S f O*  C **  - x-  ) - z^2  <A  - *\  ) ) 

° c L 

- *-  ) 

+ S\l(  T+  - n TV  ) - ^ C T*1  v >7 ov  ) 


(46) 


* S*‘CT»  - Cr-,.)1)  * i*,l(  t;  - Oz / J } <ie  dbt . 


Summand  two  tells  us  that  either  -/?=  o or  vj-  = , on  o 

Suppose  ‘/<va,  then  x**x^4=i,2.  Summand  one  then  gives  us  that 
z z 

s Xj  on  . Thus  we  have 

^ C X - 5 ) =2  > c.  . 

as  required,  and  *e  is  defined  as  the  subset  of  'c  where  x;1-*1  . 

The  last  four  summands  give  the  momentum  balance  conditions,  as 
usual,  and,  in  addition,  the  last  two  summands  imply  that  the  normal 
tractions  are  compressive  (since  ( o ).  Thus  we  have  the 
conclusions  of  Theorem  I'  and  condition  (5). 

The  analogous  set  up  for  Theorem  II'  is  accomplished  by  setting 
<r  = © in  (45)  yielding 

- O?  )*  C *1  - **  ) (47) 


for  the  integrand  of  ‘TT’l  . With  this  Eq.  (46)  becomes 
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S S { -2-  - x^))  +- 

° c ' 

- S xj-  T;  «.  \ri  T/ 

-*■  V iKj  ( T3  - CtZ  ^ cLF  rLt 

In  this  case  we  achieve  the  conclusion  of  Theorem  II'  and  condition  (5). 

Thus  to  Theorems  I'  and  II*  we  can  delete  condition  (5),  add  the 
simplifications  manifested  in  (45)  and  (47),  and  achieve  the  conclusions 
of  Theorems  I*  and  II',  respectively,  plus  condition  (5). 
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6.  Contact  Problems  for  One,  Two  and  Three-dimensional  Bodies 

The  previous  work  needs  only  trivial  modification  to  be  made 
applicable  to  contact  problems  involving  bodies  of  different  dimensions. 
There  are  many  cases  of  considerable  interest  which  fall  into  this  cate- 
gory. For  example,  models  consisting  of  a shell  and  a plate,  or  a solid 
and  a plate,  are  useful  for  the  study  of  head  impact.  The  modifications 
necessary  are  essentially  interpretative.  An  example  illustrates  this 
assertion. 

Consider  the  frictionless  Hertzian  contact  of  a three-dimensional 

V 2. 

solid  and  a two-dimensional  plate.  Let  <K>  represent  the  solid  and  '£>  the 
plate.  In  evaluating  JI  , the  <£>*  part  is  as  before  while  the  part 
would  manifest  the  particular  plate  theory  used.  The  contact  term  Xl 
(or 7 71  ) would  be  exactly  as  before.  However  note  that  ^ (or  c ) is,  in 
this  case,  also  identifiable  with  part  of  the  two-dimensional  "volume"  of 
the  plate,  rather  than  its  boundary.  Taking  variations,  everything  is  as 
before  except  that  the  term  t (or  -(n'f'zxt  ) contributes  to 
the  transverse  momentum  equation  of  the  plate,  rather  than  to  its  boundary 
conditions.  The  interpretation  of  t:  (or  - (.**)*  ) is  thus  two-fold, 
i.e.,  U is  the  normal  component  of  the  traction  vector  with  respect  to 

, as  before,  and  it  is  also  the  equivalent  normal  "body  force"  with 
respect  to  '£>*  , manifested  by  the  interaction  with  ® 1 (Fig.  7). 

This  interpretation  is  general,  namely,  for  one  and  two-dimensional 
bodies  the  contact  force  is  an  equivalent  "body  force"  which  contributes 
to  the  momentum  equations,  rather  than  the  boundary  conditions.  With  this 
interpretation  in  mind,  the  construction  of  variational  theorems,  analogous 
to  the  ones  constructed  in  Sections  2,  4 and  5,  for  the  class  of  one,  two 
and  three-dimensional  contact  problems,  is  just  a formal  deductive  exercise 
involving  only  appropriate  definitions  for  H. 
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7.  Impact 

The  previous  sections  deal  with  spatial  aspects  of  contact  problems. 

In  this  section  we  investigate  temporal  considerations,  i.e.,  those 
phenomena  which  are  unique  to  dynamic  contact  or  impact.  To  manifest 
the  problem  encountered  in  such  situations  consider  the  following 
hypothetical  situation.  Assume  that  we  are  in  the  process  of  numerically 
solving  some  impact  problem  and  suppose  that  it  is  discovered  as  we 
monitor  the  motior  of  the  bodies  that  they  impact  somewhere  in  the  time 
interval  (t-j^).  At  time  t^  we  know  the  states  of  both  bodies  and  we 
know  that  somewhere  between  and  they  have  coalesced  over  a portion 
of  their  boundaries.  Assume  for  the  moment  we  know  the  geometry  of  the 
contact  surface  fc.  The  question  which  arises  then  is  what  is  the  state 
of  k at  time  t2.  i.e.,  what  are  the  velocity  and  traction  vectors  on  k? 

It  is  necessary  to  know  this  information  to  carry  forth  the  step  forward 
time  integration.  The  quertion  though  seems  improperly  posed  without 
specifying  considerable  data  about  the  nature  of  the  impact.  To  get  a 
handle  on  things,  we  will  initially  formulate  a simple  one-dimensional 
problem  involving  the  impact  of  two  elastic  rods.  Although  this  problem 
is  trivial,  it  provides  considerable  insight  into  the  general  nature 
of  impact  of  continuum  bodies.  Since  we  are  interested  in  the  state  of 
t (in  this  case  a noint)  iumediately  after  impact,  whether  the  rods  are 
finite  or  semi-infinite  is  immaterial. 

Assume  that  the  pre-impact  states  of  the  two  bodies  are  given  by 
the  following  data: 

43  ***  , ( , f-C  > <*-'.2  (48) 
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At  impact  the  rods  coalesce  at  e,  and  for  some  finite  time  interval 
thereafter  (at  least)  xe  t is  persistent.  At  the  moment  of  impact 
shock  waves  begin  to  propagate  in  each  body.  The  space-time  picture  is 
depicted  in  Fig.  8.  As  discussed  in  section  1,  since  fc.  is  material  and 
x is  persistent,  we  have 


,f  ‘ * 

- V . - V. 


(49)* ** 


for  the  post-impact  state  (t+).  In  addition  to  (49),  the  well  known 
shock  conditions  must  hold  across  the  wave  fronts: 


C^'J  v u-  Qt  *x /<?x  )“<  3 « ( j > 


(50) 


where  l * is  the  material  velocity  of  the  shock  in  , and  [ ] is 

the  wave-front  jump  operator  which  assigns  to  a function  the  difference 
in  its  values  behind  and  in  front  of  the  wave,  i.e.,  f.  f(X,t)  - f(X  ,t) 

where  X is  a material  point  denoting  the  location  of  the  wave-front.  As 
can  be  deduced  from  Fig.  8,  the  states  into  which  the  shocks  initially 
propagate  are  the  pre-impact  states  given  by  (48),  and  the  state  at  fe, 
immediately  after  the  shocks  pass,  is  given  by  the  post-impact  state  (49). 
These  observations  in  conjunction  with  (50)  yield, 

- n/  + u- (c dy/dx £ 

PTLT  c-r  - + P'.‘ 


- (ax/px  } = c , 

- P r_ 


(51) 


*For  convenience  we  choose  the  initial  state  to  be  the  pre-impact  state, 
thus  we  need  not  distinguish  between  Cauchy  and  Piola  tractions. 

**A  consistency  condition  for  these  equations  is  that  X.'  -x?  >o. 
Otherwise  the  impact  would  not  occur. 
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The  four  Eqs.  (51)  and  constitutive  equations  relating  p"*  to  (dx/dx  T* 
yield  a formally  deterministic  system  of  six  equations  in  the  six 
unknowns  P,  U , (ax/dX)  . Thus  we  see  that  the  desired 
quantities  v and  P depend  on  the  pre- impact  states  and  material  properties 
of  both  <23  . The  precise  form  of  this  relationship  depends  upon  the 
constitutive  equations  of  the  bodies.  As  a simple  example,  assume  we 
have  linear  constitutive  equations  P“  = E~{C5x/©xT  - l } , 

E"  constant,  and  let  the  pre-impact  state  be  given  by 

-r  = V*  , 

C2>*/2>X'T  = 1 > 

P “ = O 

These  conditions,  when  inserted  in  Eqs.  (51),  lead  to: 

>✓  - 

- pju* 

P = V*  - V1  (53) 

(U*)!=  e-  /e.*  • 

Note  that  the  denominators  in  Eqs.  (53)^  ^ Present  no  problems  since 
^>7  > o , E**  > o and  ^-1  - agn  U*"  = - agn  U*  . 

This  result  is  also  appropriate  whenever  the  intensity  of  the  impact  is 
small  enough  such  that  the  non-linear  constitutive  equation  can  be 
replaced  by  its  linear  approximation  about  the  pre-impact  state.  In  this 
case  E * is  a tangent  modulus  evaluated  at  the  pre-impact  strain 

to*  /dv)_-  - l } = o .To  further  simplify,  consider  the 

case  when  both  rods  have  identical  properties  ( i . e . , pf1  , 

E = T “ , » i , 2.  ).  Then 
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P = p.U  (V‘-  V')/Z 
(uf-  e /e.  . 


(54) 


In  £qs . (54)  U is  positive,  and  since  consistency  requires  V‘-V*>  o , 
P is  compressive. 

Thus  for  the  one-dimensional  case  at  least  the  problem  of  computing 
the  post-impact  state  is  easily  achieved.  The  solution  of  (51)  for  the 
fully  non-1 inear  case  can  be  automated  as  part  of  a numerical  algorithm. 
Although  this  problem  is  trivial,  it  serves  to  indicate  that  the  post- 
impact problem,  the  solution  of  which  is  essential  in  a numerical 
algo  -ithm,  is  one  of  wave  propagation. 

In  the  analysis  of  higher  dimensional  bodies  the  solution  of  the 
post- impact  problem  becomes  greatly  complicated  due  to  the  geometric 
variety  of  impact  conditions.  However,  considerable  simplifications 
can  be  taken  advantage  of  if  one  keeps  in  m’nd  the  nature  of  the  discrete 
problem.  For  instance,  if  a certain  portion  of  the  boundaries  of  two 
bodies  have  coalesced  in  e,  each  interior  point  of  e,  at  which  the 
tangent  plane  is  well  defined,  may  bo  treated,  to  the  first  order,  as  a 
point  on  the  mating  surface  of  two  impacting  half-spaces.  As  long  as 
time  steps  are  kept  small  enough,  the  local  behavior  is  well  represented. 
The  post-impact  problem  for  the  general  case,  analogous  to  (51),  can  be 
automated  as  part  of  the  numerical  algorithm,  and  for  many  simple  cases 
can  be  solved  explicitly. 
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With  these  notions  in  mind,  let  us  return  to  the  case  of  main  interest 
in  this  report,  namely  three-dimensional  continuum  bodies.  We  shall  con- 
sider only  the  case  of  a frictionless  contact  surface  (Case  II),  and  leave 
the  solution  of  the  post-impact  problem  for  the  no-slip  case  (Case  I), 
which  is  more  difficult,  for  future  work.  With  the  proper  interpretations, 
the  one-dimensional  rod  formulation  (Eqs.  (48-54))  suffices  to  completely 
characterize  this  case.  This  is  so  because  no  tangential  motions  or 
stresses  may  be  communicated  across  a frictionless  surface,  and  thus  we 
need  only  consider  the  configuration  of  normal  incidence.  In  this  case 
the  requisite  constitutive  functionals  in  (51)  would  be  those  relating 
P04  , the  normal  Piola  stress,  to  the  normal  component  of  strain,  holding 
all  other  components  of  strain  fixed  at  the  pre-impact  values.  For 
example,  in  the  linear  isotropic  case,  E * (Young's  modulus)  in  Eqs.  (53,54) 
would  be  replaced  by  are  the  Lame  and  shear 

moduli,  respectively)  and  the  propagation  velocity  would  be  that  of 
dilatational  waves. 
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PART  II 

A NUMERICAL  SCHEME  FOR  ANALYSIS  OF 
CONTACT1- IMPAct  PROBLEMS 

8.  Numerical  Solution  of  Contact- Impact  Problems 

In  performing  numerical  computations  based  on  the  above  described 
variational  formulation  for  contact-impact  problems  we  have  employed 
three  distinct  levels  of  approximation:  (1)  a spatial  discretization 
of  the  bodies  and  contact  surfaces,  (2)  a temporal  discretization  to 
determine  the  response  of  the  discretized  bodies,  and  (3)  a numerical 
solution  for  the  resulting  system  of  nonlinear  algebraic  equations. 

In  the  following  sections  we  shall  restrict  our  attention  to  the 
Hertzian  contact  problem  described  in  Section  5.  Significant  numerical 
difficulties  are  encountered  in  the  solution  of  impact  problems;  to 
complicate  the  problem  further  by  introducing  the  additional  steps 
necessary  to  determine  the  contact  surface  maps  for  the  full  kinematically 
nonlinear  case  is  left  for  a future  study.  While  this  is  a simple  impact 
problem  in  terms  of  determining  the  contact  surface  and  the  full  power  of 
the  preceding  theory  is  neither  necessary  nor  exploited  in  its  solution, 
many  of  the  features  of  the  general  problem  are  employed  here. 
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9.  Spatial  Discretization  of  the  Bodies  and  Contact  Surface 
1 2 

The  bodies  and<©  are  discretized  using  standard  finite  element 

methods,  (e.g.t  see  [5]).  In  order  to  facilitate  the  computation  of 

i discrete  Hertzian  contact  surface  the  nodes  of  are  arranged  so 

2 

that  they  align  with  the  nodes  of  # . This  is  consistent  with  the 
notions  of  condition  3 of  Section  5 and  ensures  that  during  determination 
of  the  approximation  to  the  contact  surface  contiguous  nodes  of  the  two 
bodies  will  meet.  Thus,  the  simulation  of  the  contact  surface  is  trivial. 
The  development  of  a numerical  model  for  Hertzian  contact  problems  is 
based  upon  the  form  of  Theorem  II 'which  uses  (47)  for  the  integrand  of 
7T7.  For  numerical  computations  we  introduce  the  displacement  vector 
u such  that 

k = X”  + u (55) 

For  a compatible  finite  element  displacement  field  the  integrand  of 77? can 
be  approximated  by  taking  as  the  product  of  CO  and  S(z-Zi) 

(i.e.,  Dirac  delta  functions  in  space).  This  corresponds  to  taking  z 1 
as  "concentrated  nodal  loads"  which  are  the  generalized  forces  of  the 
contact  pressure.  With  this  discretization  we  can  describe  pseudo 
contact  elements  between  each  pair  of  candidate  contact  nodes.  Let 
these  nodes  be  denoted  as  ( f and  the  generalized  force  as 

k 

then 


= \ - u,l.w  * *:  - >t  (56) 

o 

where  (i)  are  the  set  of  candidate  contact  nodes  which  span  c ; are 
the  nodal  displacements  in  the  z3  direction  and  are  the  nodal  coordinates 
of  the  candidate  contact  nodes.* 

*We  assume  here  that  3 is  the  direction  nominally  normal  to  the  contact 
surfaces,  e.g.,  see  Fig.  5. 
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Use  of  the  finite  element  method  in  Theorem  Il'with'TTt,  given  by 
(56)  produces  a set  of  nonlinear  second  order  ordinary  differential 
equations  which  together  with  the  impenetrability  conditions  define 
the  discretized  contact  impact  problem.  These  equations  take  the  form: 

Mu  + K Cu)  = R , (57) 

where  M is  the  usual  finite  element  mass  matrix, K represents  the  elastic 
stiffness  forces  together  with  the  contact  terms, R is  the  set  of  general- 
ized forces  resulting  from  boundary  tractions  and  u is  the  set  of  time 

2 

dependent  nodal  displacements  (which  also  include  the  (£,)).  For  inelastic 
materials  Theorem  Il’can  be  extended  by  treating  the  first  variation  as  a 
Galerkin  method  (principle  of  virtual  work)  and  replacing  the  elastic 
constitutive  model  by  more  general  theories,  e.g.,  viscoelastic,  elasto- 
plastic,  viscoplastic,  etc.  In  this  case 

— K(u,u)  (58) 


in  (57). 
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10.  Temporal  Discretization 

A temporal  discretization  of  the  second  order  ordinary  differential 
equations  which  result  from  a finite  element  spatial  discretization  of 
the  contact-impact  problem  is  accomplished  herein  by  using  the  Newmark 
family  of  methods  [6].  The  Newmark  family  of  methods  is  a one-step 
integration  method  with  two  free  parameters  which  can  be  used  to  control 
stability  and  numerical  damping.  The  method  is  essentially  a difference 
method  in  time.  The  behavior  of  the  method  for  linear  elasto-dynamics 
problems  is  discussed  in  [6,7].  The  algorithm  is  given  by 

u =■  V-* ^ ~ - (3  )AtV„  + p 6tl  Ci  , (59 

and  ^ u.„  * ir  at  ui 


where  u „ = u.  (t„)  y at  - , and  /3  , f 

are  the  two  parameters.  For  linear  problems  if=.5*s  = .5  produces 
no  artificial  viscosity  and  /3  * ^ C'4  * )a  produces  unconditional 
stability  (i.e.,  the  method  is  stiffly  stable).  Such  generalization 
is  not  possible  for  nonlinear  problems  and  during  solution  it  may  be 
necessary  to  monitor  the  solution  for  any  signs  of  instability.  In  (59) 
/3  =*  o . produces  an  explicit  method  for  ii and  if  hi  is  diagonal 
(lumped  mass)  with  K and  independent  of  u the  solution  can  be  advanced 
without  solving  a large  set  of  simultaneous  equations;  for  all  other 
cases  the  method  is  implicit  and  equations  must  be  solved.  Solution  of 
(59)  1 for  um,  in  terms  of  the  solution  at  tn  and  un>)gives 


^ = -l-  m 

~ ft  at 


- x 

/Sat 


(60) 
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which  can  also  be  used  in  (59)2  to  express  the  velocity  in  terms  of  the 
solution  at  t„  and  u.^(.  Since  in  this  process  we  divide  by  /3  and  At 
it  is  no  longer  possible  to  consider  zero^  or  zero  time  steps. 


11 . Solution  of  the  Nonlinear  Algebraic  Problem 

Use  of  the  Newmark  method  in  (57)  (including  (58))  yields  the  set 
of  nonlinear  algebraic  equations: 

-*•  *5  ,y-„  ) = R +-M  A„  > (61) 


where 


A 


_J 

/3^C^ 


/B  At 


r ^ ^ - 

^ a./3  -7 


A Newton-Raphson  iterative  solution  to  this  set  of  equations  can  formally 
be  constructed,  giving: 


( !__  M * 


- m 
14 


(62) 


wnere  is  the  effect  of  loads  varying  with  the  deformation 

and 


C 2>-K  )L>  = 


(63) 


is  the  tangent  stiffness  matrix.  The  coefficient  to  aJ"\'s  generally 
called  the  Jacobian  matrix  of  the  Newton-Raphson  iteration.  The 
solution  is  advanced  by  taking 


<»-!>  (•*■>  t*> 

V A U.  % 


(64) 


and  iterating  until  a norm  of  the  solution  satisfies 
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where  e is  some  small  positive  error  tolerance.  In  the  work  reported 
here  the  norm  II  II  is  taken  as  tta  Euclidian  norm 

II  * (I  = **  yA  , (66) 

and  the  load  vector  R is  assumed  to  be  independent  of  a.  For  stable 
elastic  materials  the  resulting  tangent  stiffness  is  then  symmetric  and 
positive  definite,  consequently,  standard  direct  solution  methods 
normally  employed  in  the  solution  of  linear  finite  element  problems  can 
be  used.  For  inelastic  materials  or  deformation  dependent  loads  the 
tangent  stiffness  may  be  asymmetric.  In  these  cases  some  special  methods 
may  be  necessary  to  effect  a solution. 
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12.  Discretized  Impact  Conditions 

In  the  previous  numerical  development  o'  has  been  defined  by  discrete 
points  which  correspond  to  nodes  along  the  boundaries  of  and  <6  . 

When,  during  the  course  of  advancing  the  solution  in  time,  any  one  of 
these  points  violates  the  impenetrability  condition  a re-solution  must 
be  obtained  in  which  the(t^  are  now  non-zero  and  the  u.*  satisfy  the 
impenetrability  condition.  Some  control  and  monitoring  are  required  to 
effect  this  in  a computer  program.  In  addition  to  satisfying  these 
conditions,  the  impact  relations  denoted  in  Section  7 must  be  invoked. 

In  the  present  study  these  conditions  are  applied  to  the  solution  at  the 
end  of  a time  step  in  which  points  first  go  into  contact.  Accordingly 
we  compute  from  (50)* 


(67) 


1 2 

and  assign  this  value  to  the  appropriate  node  of  </3  andtB  . 

To  determine  the  solution  vector  a at  tn4.,  we  have  solved  the  set 
of  equations  (61).  As  described  above  the  shock  conditions  are  then  used 
to  determine  the  value  of  the  velocity  at  time  tn+1  for  all  points  which 
have  come  into  contact  during  the  time  interval.  In  order  to  get  a con- 
sistent solution  at  these  points  we  must  modify  the  accelerations  and 

contact  force  to  reflect  the  shock  conditions.  This  is  accomplished  by 

1 2 

re-solving  the  equilibrium  conditions  of  <9  and  <-9  at  point  i.  The 
expanded  forms  of  the  appropriate  equations  are: 


*The  ( )_  denotes  a value  which  is  computed  before  impact,  whereas  ( )+ 
denotes’the  value  after  impact. 
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Mlu_  v Kl(u)  - (£i)l  = R‘  , 

♦ k'(u)  = R2  . 


(68) 


For  nodes  which  have  come  into  contact  we  must  enforce  the  condition  on 
acceleration 


ul 


(69) 


2 

and  compute  the  contact  force  (£^)  +.  The  solution  for  these  is  obtained 
from 

M'u.  * K‘Cy)  ♦ (&l)\  - , 

and 

> kg*)  * Ce-\  = • 

These  are  two  equations  in  two  unknowns  which  can  be  solved  for  the  Ci+ 
and  (t^)£.  If  is  independent  of  velocity  the  stiffness  forces 

and  R^will  remain  unchanged  during  the  impact,  hence  we  can  solve  the 
simpler  problem 

M'u+  - (Ll£  = Kl'C^  - 

Mzu+  - COt  = ~ C tS_ 

whose  solution  is 

ul  + « M u_  M.  UL  _ 

M 1 * M1  ’ 


*C£t)t  * *COt  + m‘  (ul  - ) 

- M!  (u  - u,  ) . 


and 
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I 


i 


I 


This  completes  the  numerical  specification  of  the  solution  at  tn+1 ; this 
solution  process  is  now  repeated  for  each  of  the  succeeding  time  steps. 

At  this  point  it  is  important  to  compare  the  solution  procedure  for 
impact  of  a continuum  discretized  by  a finite  element  method  with  the 
solution  procedure  for  a physically  discrete  body,  i.e.,  a body  composed 
of  mass  points  joined  by  massless  elastic  springs.  Both  problems  may  be 
described  by  algebraic  equations  of  the  form  of  (57).  The  impenetrability 
condition  is  also  identical.  The  impact  conditions,  however,  are  different. 
For  the  discretized  continuum  the  procedure  is  described  above.  The  study 
of  the  impact  of  mass  points  is  considered  in  elementary  mechanics  books, 
e.g.  [6],  The  impact  of  two  mass  points  is  described  by  impulsive 

motion  such  that  at  t the  velocities  of  the  two  mass  points  are  V1  and 

2 1 2 
V ; after  impact  at  time  t+,  the  two  points  have  velocities  V+  and  V+  . 

1 2 

The  two  points  will  not  in  general  stay  in  contact  (i.e.,  V+  ^ V+  ) 

1 2 

but  will  rebound.  The  conditions  used  to  compute  the  V+  and  V+  are: 
Balance  of  Momentum* 

t-V  { V'}  M1  ■(  v1}  = o , (7i ) 


and  use  of  an  equation  involving  the  "coefficient  of  restitution",  e : 


V*  - 


(72) 


For  e=t  energy  is  conserved  whereas  for  e=o  the  points  "stick"  and 
energy  is  dissipated.  We  must  conment  in  passing  that  (72)  is  the  energy 


* 


-C=fco>  = fct.)  - f(t.) 
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equation  in  disguise.  To  see  this  we  can  write  the  jump  conditions  for 
energy  as 


+ j.  mU(vi)‘>  = -o> 
z z 


(73) 


The  term  {J  > can  exist  only  if  other  energies  are  dissipated  during 
the  jump.  We  rewrite  (73)  by  using 

- [v]<v>  , 

where 

<V>-  i(N4*  W.*)  . (74) 

Use  of  the  momentum  equation  (71)  then  gives,  after  dividing  by 


<V'>  - <V->  - -[■>> 

M*  IV'} 

or  after  recollecting  terms  and  dividing  by  - V2  ) we  obtain: 

~ = 1 - -Q  > . (75) 

v* _ v2  - v.1) 

The  significance  of  the  coefficient  of  restitution  then  is  associated  with 
the  right  hand  side  of  (75). 

It  is  clear  from  the  above  developments  that  the  numerical  simulation 
of  the  discretized  continuum  and  the  physically  discrete  system  involve 
two  distinct  methods  for  treating  the  impact  conditions.  It  is  imperative 
then  to  associate  the  correct  method  for  the  problem  at  hand.  In  the 
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present  study  we  are  interested  in  the  impact  of  continua,  and  in  this 
case  we  shall  employ  the  discrete  shock  condition  to  effect  the  solution. 
This  a priori  assumes  that  the  response  we  are  computing  involves  a 
time  scale  associated  with  wave  propagation  problems.  Consequently,  we 
cannot  expect  the  computation  procedure  for  advancing  the  solution  in 
time  to  be  accurate  if  we  take  time  steps  greatly  in  excess  of  transit 
times  through  each  body.  In  this  context  it  may  be  important  to  consider 
an  "explicit"  time  integration  procedure  in  future  work.  The  stability 
restrictions  may  be  too  severe  to  make  this  feasible. 
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PART  III 

FEAP  74  - A COMPUTER  PROGRAM  FOR 

SOLUTION  Of  CONTACT- iHpACT  pBMeEs 

13.  Development  of  a Contact-Impact  Model  for  FEAP 

In  order  to  incorporate  an  ability  to  compute  solutions  to  contact- 
impact  problems  using  a finite  element  method  as  described  above  it  is 
necessary  to  have  available  a computer  program  which  can  solve  the 
nonlinear  equations  of  motion  given  by  (61).  The  computer  program  FEAP 
is  a general  program  to  solve  finite  element  problems.  The  program  has 
a capability  of  solving  both  quasistatic  and  dynamic  oroblems  and  can 
incorporate  several  types  of  elements  simultaneously.  The  nonlinear 
capabilities  required  for  the  solution  of  contact- impact  problems  have 
been  incorporated  into  FEAP  and  currently  includes  the  user  options 
(see  Input  Instructions): 

(1)  Selection  of  quasistatic  or  dynamic  option:  The  dynamic  option  will 
integrate  the  equations  of  motion  using  the  one-step  Newnark  method 

to  advance  the  solution  in  time.  Quasistatic  analysis  is  accommodated 
by  any  one-step  algorithm.  The  algorithm  employed  is  incorporated 
into  each  element  routine  and  thus  id  defined  by  the  developer  of 
each  element.  Impact  problems  require  description  of  the  contact 
surface  and  wave  speeds. 

(2)  Selection  of  the  nonlinear  method  to  advance  the  solution:  Options 

include: 

(a)  No  iterations  in  each  time  step.  Unbalanced  forces  at  each 
time  are  added  to  the  next  time  step. 

(b)  Iterations  in  each  time  step  to  achieve  a balance  of  force  within 
each  time  step.  In  this  option  the  user  can  select  to  reform 
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the  Jacobian  matrix  for  each  iteration  or  only  at  the  first 
iteration  in  each  time  step. 

In  the  impact  problems  solved  to  date  it  has  been  necessary  to  use 
the  general  form  of  the  Newton-Raphson  algorithm.  This  includes  a complete 
formating  and  factoring  of  the  Jacobian  matrix  for  each  iteration  of  each 
time  step  in  the  analysis.  If  the  method  dercribed  herein  is  to  become 
computationally  effective  improvements  in  the  computer  program  are  para- 
mount. Undoubtedly  the  most  important  aspect  in  reducing  computer  times 
is  to  introduce  a substructuring  system  so  that  the  highly  nonlinear 
equations  in  the  vicinity  of  the  contact  surface  can  be  isolated  from  the 
remainder  of  the  bodies.  This  will  normally  involve  only  a small  number 
of  equations  in  the  total  system  of  (62).  The  solution  of  a large  finite 
element  problem  will  generally  concentrate  the  computer  solution  time  in 
the  forming  and  factoring  of  the  tangent  stiffness  matrix.  The  fewer 
times  that  it  is  necessary  to  perform  this  costly  step  the  more  efficient 
tne  solution  algorithm.  Substructuring  can  be  used  then  to  restrict  the 
part  of  the  equations  which  must  be  formed  and  factored  often,  and  thus 
greatly  reduce  the  computer  costs  in  analyzing  impact  problems. 

The  version  of  FEAP  which  can  currently  be  used  to  analyze  contact- 
impact  problems  includes,  in  addition  to  the  nonlinear  Newton-Raphson 
iterative  algorithm,  a new  special  contact-impact  element  and  a new  sub- 
routine to  describe  impact  surfaces  and  the  discrete  shock  conditions 
described  in  Section  12.  These  are  described  in  the  following  sections. 
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14.  Contact  Element  for  Hertzian  Contact 

The  contact- impact  element  which  has  been  developed  is  called  ELMT05 

and  can  be  used  along  any  coordinate  direction.  As  developed  it  cannot 

be  used  a long  normals  which  are  in  non-coordinate  directions.  The 

development  of  the  contact  element  assumes  that  within  the  framework  of 

1 2 

linear  elasticity  theory  a node  on  .)  will  impact  on  a node  off'3  . In 

2 

using  this  contact  element  we  shall  assume  that  the  contact  surface  on  ' 
is  located  at  larger  coordinate  values  than  the  contact  surface  of  -J1 . 

The  contact  elenent  is  described  by  three  nodes.  Node  1 is  associated 
with.d.  Node  3 is  asso  ated  with«£,  and  Node  2 is  used  as  storage  for 
the  contact  force  (6)  . The  user  can  select  the  direction  of  contact 
motion  by  specifying  the  degree  of  freedom  of  the  nodal  unknowns  to 
which  the  contact  is  to  be  measured;  this  must  agree  with  the  physical 
direction  of  the  element  (see  Fig. 9 ).  The  degree  of  freedom  for  the 
contact  element  is  specified  during  the  MATERIAL  data  input  and  consists 
of  a single  card  in  15  format.  The  element  nodes  are  described  along 
with  all  other  elements  according  to  Section  4 of  the  Input  Instructions. 
The  Node  order  as  shown  in  Fig.  9 must  be  observed. 
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15.  Impact  Surface  Description 

The  definition  of  the  Impact  surface  includes  a list  of  all  elements 
on  the  contact  surface  together  with  the  degree  of  freedom  describing  the 
direction  of  contact  motion  (as  described  above).  In  addition,  the  pro- 
duct of  mass  density  and  wave  velocity  (always  a positive  number)  for  each 
body  is  input.  This  assumes,  currently,  that  (1)  each  contact  surface 
belongs  to  a linear  material,  and  (2)  the  same  material  exists  along  all 
of  the  contact  surface.  This  data  need  be  prescribed  only  for  impact 
problems,  quasistatic  contact  problems  do  not  require  this  data  since 
no  velocity  or  acceleration  computations  are  performed  in  this  class  of 
problems.  Data  to  be  input  for  the  impact  surface  is  given  in  Table  I. 


Table  I - Impact  Surface  Data 


CARD 

1) 

(6X,A6) 

COL. 

7 to  12 

CARD 

2) 

(2F10.0) 

COL. 

1 to  10 

pU 

COL. 

11  to  20 

?U 

CARD 

3) 

(15) 

COL. 

1 to  5 

CARD 

4) 

(215) 

Repeat  NLIST  times 

COL. 

1 to  5 

Must  contain  CONTAC 

of  body  1 
of  body  2 

NIIST,  number  of  elements  on  contact 
surface 


COL.  6 to  10 


Contact  element  number 

Degree  of  freedom  of  this  contact  element 
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16.  Example  Problems 

Two  example  problems  are  included  to  illustrate  the  characteristics 
of  the  methodology  and  the  associated  computer  program  described  above 
for  Hertzian  contact  problems.  The  first  problem  is  a quasistatic  contact 
problem  which  is  used  to  demonstrate  the  ability  of  the  computer  program 
to  compute  an  evolving  contact  surface.  The  second  problem  will  demon- 
strate the  ability  of  the  program  to  properly  model  the  temporal  response 
of  an  impact  problem. 

To  model  a problem  in  which  a contact  surface  will  change  under 
different  load  levels  we  consider  two  beams  with  an  initial  parabolic 
curvature.  A symmetric  configuration  is  analyzed  and  the  resulting  finite 
element  model  is  shown  in  Fig.  10.  Each  element  is  nominally  one  unit  by 
uiie  unit.  The  gap  at  the  load  end  is  initially  0.5  units.  The  material 
properties  used  are  E = 500  and  v = 0.  The  load  P is  applied  as  shown 
and  allowed  to  increase  linearly  in  time.  The  problem  then  is  to  deter- 
mine the  contact  surface  at  various  load  levels.  In  order  to  eliminate 
a singularity  in  the  system  of  equations  it  was  necessary  to  permanently 
attach  the  two  nodes  at  the  symmetry  axis  of  the  contact  surface.  All 
other  nodes  along  the  boundaries  between  the  two  bodies  are  assumed  to 
be  possible  cont'ct  points  and  contact  elements  are  assigned  between 
vertical  nodal  pairs.  The  load  was  varied  from  0.2  to  0.8  in  increments 
of  0.1  and  the  computed  contact  surface  and  forces  were  computed.  These 
are  given  in  Table  II.  The  deformed  shape  at  a load  of  0.4  is  also 
shown  in  Fig.  10  as  a dotted  form.  The  attached  node  at  the  center  has 
influenced  the  solution  at  loads  above  0.3  since  the  contact  pressure 
there  is  tensile  (neqative).  The  force  is  small  and  should  not  greatly 
affect  the  actual  contact  region  computed.  As  the  load  increases  the 
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contact  surface  moves  toward  the  load.  This  is  conceptually  correct  since 
if  the  beams  were  modeled  according  to  Euler-Bernoull  i theory  the  contact 
force  would  be  a point  load  which  gradually  moves  from  the  center  to  the 
outet  edge  according  to  the  relationship  (using  the  above  values  for  sizes 
and  material  properties) 


X 


= 5C'- 


This  relation  predicts  that  the  contact  point  will  be  non-zero  only  after 
P exceeds  1/6.  The  finite  element  model  is  in  qualitative  agreement  with 
this  beam  theory,  but  since  shear  deformations  are  included  the  finite 
element  solution  gives  a distributed  load  on  the  contact  surface.  It  is 
interesting  to  also  note  that  the  contact  force  over  the  center  of  the 
beams  is  zero, just  as  in  the  beam  theory. 

Table  II  - Contact  Forces 


LOAD 

X-COORDINATE 

BEAM 

THEORY -X 

0 

^1 

2 

3 

4 

5 

0.2 

0.2 

- 

- 

- 

- 

- 

0.83 

0.3 

.07 

.23 

- 

- 

- 

- 

2.22 

0.4 

-.01 

.07 

.34 

- 

- 

- 

2.92 

0.5 

-0.00 

.02 

.23 

25 

- 

- 

3.33 

0.6 

-.01 

- 

.09 

51 

.01 

- 

3.61 

0.7 

-.01 

- 

.07 

45 

.19 

- 

3.81 

0.8 

-.01 

- 

.05 

38 

.38 

- 

3.96 

This  problem  demonstrates  that  the  computer  program  can  model  the 
evolution  of  a contact  surface.  Of  particular  importance  is  to  note 
that  as  the  load  increases  the  program  can  both  attach  and  detach  a 
contact  point.  This  is  an  essential  requirement  for  the  analysis  of 
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impact  problems  as  is  shown  in  the  next  problem. 

As  a simple  example  we  consider  the  impact  against  a rigid  wall  of 
a finite,  linear  elastic  rod  traveling  at  constant  velocity.  The  rod  has 
a modulus  of  elasticity  E of  100,  and  a mass  density  p of  0.1.  The  arrival 
velocity  is  taken  to  be  0.1  (units  may  be  assigned  in  any  convenient  system). 
The  rod  is  taken  to  be  10  units  long  and  is  divided  into  10  elements  plus 
one  contact  element  as  shown  in  Fig.  11.  At  time  zero  the  rod  is  just 
arriving  at  the  wall.  The  exact  solution  predicts  a contact  duration  of 
0.2  time  units.  This  corresponds  to  the  time  required  for  a wave  to  travel 
from  the  contact  point  to  the  left  end  and  back  to  the  contact  point  at 
which  time  the  rod  will  part  from  the  wall.  The  problem  was  analyzed 
using  FEAP  with  time  steps  of  0.01  unit  (transit  time  across  an  e^ent) 
and  the  rod  remains  in  contact  until  time  0.20  units  and  has  rebounded  at 
time  0.21.  Thus  the  program  can  predict  accurately  the  contact  duration 
of  the  rod.  The  finite  element  solution  obtained  is  compared  with  the 
exact  solution  in  Fig.  12.  The  agreement  of  stresses  and  contact  force 
is  good.  The  largest  discrepancy  exists  in  defining  the  shock  front, 
which  is  "smeared"  by  the  finite  element  method  and  ordinary  differential 
equation  solution  method  used  here.  This  is  the  same  type  of  solutions 
which  are  commonly  obtained  with  numerical  solutions  of  this  type  even 
without  impact.  Solutions  such  as  the  impact  shocks  generated  are  probably 
one  of  the  most  difficult  responses  to  accurately  calculate  by  a finite 
element  method. 
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Figure  11 
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17.  Closure  and  Recommendations  for  Future  Work 

In  the  preceding  sections  we  have  presented  a theory  for  contact- 
impact  problems  together  with  the  numerical  development  of  a Hertzian 
contact- Impact  model.  The  computer  program  FEAP  74  has  been  modified 
to  include  the  model  and  has  successfully  solved  a contact  problem  and 
an  impact  problem.  The  work  reported  herein  must  be  considered  initiatory; 
the  general  theories  and  their  numerical  implementation  have  not  been 
completed.  The  problem  is  of  such  a complicated  nature  and  the  literature 
existing  prior  to  this  study  was  so  meager  that  we  consider  it  fitting 
to  document  the  work  completed  thus  far. 

We  have  attempted  to  qualify  each  stage  of  the  development  throughout 
the  report,  however.  It  may  be  fitting  to  reiterate  future  work  which  we 
consider  to  be  essential  for  numerical  models  to  be  effective  and  efficient 
tools  for  predictive  analyses. 

(1)  The  restriction  of  Hertzian  type  contact  must  be  removed.  This 
Involves  the  non-trivial  task  of  finding  appropriate  numerical 
methods  to  handle  the  maps. 

(2)  Improved  methods  for  solving  the  set  of  nonlinear  algebraic 
equations  must  be  found.  We  have  suggested  two  methods  which 
should  be  considered:  (a)  Substructure  the  problem  about  the 
contact  regime  so  that  a more  efficient  forming  and  factoring 
of  the  tangent  stiffness  can  be  performed;  and  (b)  Since  the 
impact  problem  is  a wave  propagation  problem  an  explicit 

time  integration  of  the  equations  of  motion  should  be  explored. 

In  complex  situations  the  explicit  integration  method  may  have 
severe  stability  limitations  which  could  make  it  unacceptable. 
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(3)  Methods  of  utilizing  the  shock  conditions  need  to  be  explored 
further.  We  have  noted  some  peculiar  anomalies  when  the  bodies 
separate.  These  appear  to  be  caused  by  a shock  like  separation 
phenomena. 

(4)  When  the  wave  propagation  property  of  the  impact  problem  is 
ignored  by  taking  time  steps  greatly  in  excess  of  the  transit 
times  in  a body  the  computed  response  is  meaningless.  Under 
such  situations  the  bodies  rebound  within  a single  time  step. 
Currently  the  rebound  velocity  is  much  too  large.  When  the 
shock  conditions  are  used  for  a class  of  problems  where  the 
response  desired  is  in  the  target  instead  of  in  the  impactor, 
it  may  be  expedient  to  take  a large  time  step.  Methods  should 
be  explored  to  accomplish  this  capability. 

The  above  recommendations  for  future  work  should  in  no  way  minimize 
what,  has  been  accomplished  by  the  present  study.  For  the  first  time  a 
contact- impact  theory  in  the  form  of  a variational  problem  has  been  pre- 
sented in  a general  form.  This  formulation  was  motivated  by  the  fact  that 
numerical  solutions  would  be  obtained  by  a finite  element  method.  In 
addition  the  necessary  foundation  for  the  numerical  solution  has  been 
thought  out  and  within  this  context  a computer  program  has  been  developed 
for  Hertzian  contact- impact  problems. 

The  implementations  considered  here  have  produced  results  which  are 
hopeful  signs  for  the  eventual  success  of  the  more  general  impact  problems. 
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APPENDICES 


A.  Input  Instructions  for  Contact/Impact  Problems 

In  order  to  analyze  contact/impact  problems  in  FEAP,  users  must 
prepare  the  data  for  a time  dependent  analysis.  This  will  include 
the  following  Data  Type  Identification  Cards  (see  Section  1, 
Appendix  B): 

FEAP  74 
MATERIAL 
NODAL 
ELEMENT 

CONTACT  (for  impact  problems  only) 
loadings 

INITIAL  CONDITIONS  (if  non-zero) 


VISCOE  (for  quasi-static  contact  problems) 


or 

IMPLICIT  (for  impact  problems) 

In  performing  the  necessary  solution  to  (62)  the  full  Newton-Raphson 
method  must  be  employed;  this  is  controlled  by  the  data  in  col.  76-80  of 
the  first  card  following  the  VISCOE  or  IMPLICIT  card,  and  consists  of  a 
negative  number  (negative  uses  full  Newcon-Raphson  iteration  with  the 
absolute  value  of  the  number  giving  the  number  of  iterations  to  be  per- 
formed before  going  to  the  next  time  step).  As  an  example  of  the  required 
input  data  Table  A shows  the  input  data  used  for  the  impact  problem 
reported  in  Section  16. 


. 14  OUTPUT 

FE AP<M  * * ONF  ? IfENS  IONAL  BAR  CONTACT  PPOBLt'M 
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B.  Input  Instructions  for  FEAP  74 

The  Input  instructions  for  the  description  of  a finite  element 
mesh,  together  with  the  initial  and  boundary  conditions,  is  described 
by  subroutine  MANUAL  listed  on  the  following  pages.  The  input  of  the 
contact  surface  for  Impact  problems  Is  described  in  Table  I of  this 
report. 

The  description  of  material  properties  for  the  contact  element  is 
described  in  Section  14.  For  material  properties  for  other  elements 
in  FEAP  special  input  instructions  must  be  supplied. 
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C.  Listing  of  the  Contact/ Impact  Subroutines  Added  to  FEAP  74 

The  listings  for  the  subroutines  which  are  added  to  FEAP  74  for 
the  contact/impact  theory  described  herein  are  given  below. 
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